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We propose a new method to solve the Hartree-Fock-Bogoliubov equations for weakly bound 
nuclei, which works for both spherical and axially deformed cases. In this approach, the quasi- 
particle wave functions are expanded in a complete set of analytical Poschl-Teller-Ginocchio and 
Bessel/Coulomb wave functions. Correct asymptotic properties of the quasiparticle wave func- 
tions are endowed in the proposed algorithm. Good agreement is obtained with the results of the 
Hartree-Fock-Bogoliubov calculation using box boundary condition for a set of benchmark spherical 
and deformed nuclei. 
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I. INTRODUCTION 

The study of nuclei far from stability is an increasingly 
important part of contemporary nuclear physics. This 
topic is related to newly created radioactive beams facili- 
ties, allowing more experiments on nuclei beyond the sta- 
bility line. The new experimental opportunities on nuclei 
with extreme isospin ratio and weak binding bring new 
phenomena which inevitably require a universal theoret- 
ical description of nuclear properties for all nuclei. The 
current approach to the problem is the nuclear density 
functional theory which implicitly rely on Hartree-Fock- 
Bogoliubov (HFB) theory, unique in its ability to span 
the whole nuclear chart. 

The HFB equations can be solved in coordinate space 
using box boundary condition P, Q • This approach (ab- 
breviated HFB/Box in this paper) has been used as a 
standard tool in the description of spherical nuclei 
Its implementation to systems with deformed equilibrium 
shapes is much more difficult, however. Different ap- 
proaches have been developed to deal with this problem, 
such as the two-basis method [1, |^, @, the canonical- 
basis framework 0, H, and basis-spline techniques in 
coordinate-space calculations developed for axially sym- 
metric nuclei [13, [IH ■ These algorithms are precise, but 
time-consuming. 

Configuration-space HFB diagonalization is a useful 
alternative to coordinate-space calculations whereby the 
HFB solution is expanded in a complete set of single- 
particle states. In this context, the harmonic oscillator 
(HO) basis turned out to be particularly useful. Over 
the years, many configuration-space HFB codes using the 
HO basis (abbreviated HFB/HO) have been developed. 



employing either the Skyrme or the Gogny effective in- 
teractions [l^. Tsl . [H, [H, [H, |T3] , or using a relativistic 
Lagrangian 19| in the context of the relativistic Hartree- 
Bogoliubov theory. In the absence of fast coordinate- 
space methods to obtain deformed HFB solutions, the 
configuration-space approach has proved to be a very fast 
and efficient alternative allowing large-scale calculations 

[a El- 

Close to drip lines, however, the continuum states start 
playing an increasingly important role and it becomes 
necessary to treat the interplay of both continuum and 
deformation effects in an appropriate manner. Unfor- 
tunately, none of the existing configuration-space HFB 
techniques manage to incorporate continuum effects. 

The goal of the present work is to find an efficient nu- 
merical scheme to solve HFB equations for spherical and 
axially deformed nuclei, which properly takes the contin- 
uum effects into account. We will denote this problem as 
continuum HFB (CHFB). Aiming at treating spherical 
and deformed nuclei on the same footing, we rely on the 
configuration-space HFB approach. 

The HO basis has important numerical advantages; for 
example, the use of the Gauss-Hermite quadrature al- 
lows fast evaluation of matrix elements. On the other 
hand, its Gaussian asymptotics prevents from expand- 
ing systems with large spatial extension, such as halo 
nuclear states. This problem can be successfully fixed 
by using the transformed HO basis (THO) ^20,]. The 
latter transforms the unphysical Gaussian fall-off of HO 
states into a more physical exponential decay. Neither 
HO nor THO bases, however, are able to provide proper 
discretization of the quasiparticle continuum. This has 
repercussions already at the HFB level, for which the 



2 



HO and THO bases cannot reproduce simultaneously all 
asymptotic properties of nuclear densities (see Sec. N} . 
While this shortcoming is obvious for the HO basis, it also 
arises for the THO basis because the latter can provide 
only one type of asymptotic form, i.e. the one inserted 
in the scaling function defining the THO wave functions 
[l3| ■ Hence, the THO basis fails to reproduce asymptotic 
properties, as asymptotic behavior is different for respec- 
tive channels: proton and neutron, normal and pairing 
densities, different angles for the deformed case. In fact, 
differences between calculations using the THO and the 
coordinate-space bases have been noticed in pairing prop- 
erties of nuclei (see Sec. IVland Ref. [13). This indicates 
that THO calculations may not always be fully accurate 
even in the nuclear region and necessitate careful check 
of obtained results. For the aim of carrying out quasipar- 
ticle random phase approximation (QRPA) calculations 
with the HFB quasiparticle representation, the HO and 
THO bases are very likely to be insufficient as they can- 
not provide accurate quasiparticle wave functions in the 
continuum region. 

Obviously, a more practical basis is needed. The 
Gamow Hartree-Fock (GHF) basis [l^l would be appro- 
priate, as it has been demonstrated that it can provide 
the correct asymptotic of loosely bound nuclear states. 
However, it implies the use of complex symmetric matri- 
ces. Moreover, the presence of basis states which increase 
exponentially in modulus leads to numerical divergences, 
unless the costly two-basis method is employed 

As we plan to consider bound HFB ground states only, 
it is more advantageous numerically to employ Hermitian 
completeness relations, whose radial wave functions are 
real. They are either bound, thus integrable, or oscillate 
with almost constant amplitude, so that we are free from 
the numerical cancellation problems associated with the 
Gamow states. It should be stressed that we can generate 
a Gamow quasiparticle basis using the HFB potentials 
thus obtained. We can then describe resonant excited 
states by means of the quasiparticle random phase ap- 
proximation representing the QRPA matrix elements in 
terms of the Gamow quasiparticle basis. This serves as 
an interesting subject for future investigation. 

One could expect that the employment of the spherical 
Hartree-Fock (HF) potential to generate the real contin- 
uum HF (CHF) complete basis would solve the problem. 
Unfortunately, the CHF basis is not numerically stable 
due to the presence of resonances in the vicinity of the 
real continuum. The continuum states lying close to a 
narrow resonance are rapidly changing, so that a very 
dense continuum discretization around this resonance is 
necessary to accurately represent this energy region. Im- 
portant numerical cancellations would occur as contin- 
uum wave functions become very large in amplitude close 
to narrow resonances. 

To overcome this difficulty, we adopt a technique based 
on the exactly solvable Poschl-Teller-Ginocchio (PTG) 
potential [1^. The spherical HF potential, seemingly 
the best candidate to generate a rapidly converging basis 



expansion, but providing numerically costly GHF bases 
or unstable CHF bases, is replaced by a PTG potential 
fitted to the HF potential if the latter give rise to resonant 
structure. It will be shown that the narrow resonant 
states of the GHF basis will become bound in the PTG 
basis, so that its scattering states will have no rapid phase 
shift change, a necessary condition for numerically stable 
continuum discretization. As a result, we obtain a very 
good basis for HFB calculations. We call this approach 
HFB/PTG. 

To test the feasibility of this new method, we have per- 
formed numerical calculations for spherical Ni isotopes 
near the drip line, ®^Ni - ^"^Ni, for a strongly deformed 
nucleus ^^°Zr, and two HFB solutions for ^''Mg with dif- 
ferent, prolate and oblate, deformations. Good agree- 
ment with THO calculations is obtained. 

The paper is organized as it follows. The HFB/PTG 
algorithm is described in Sec. [Hi while the method used 
to generate the PTG basis is formulated in Sec. IIIII 
Asymptotic properties of the HFB quasiparticle wave 
functions are discussed in Sec. IIVI Results of numeri- 
cal calculation are presented in Sec. [V] Brief summary 
and conclusions are given in Sec. IVII Some technical de- 
tails related to the PTG basis and calculation of matrix 
elements are collected in Appendices. 



II. THE HFB/PTG APPROACH 

Our aim is to develop an efficient method of solving 
the CHFB equation 



h{r(7, rV) - A 
h{ra, r'cr') 

/C/(i?,rV) 
H^(i?,rV) 



h{rcr, rV) 
-h{r(T, r'o-') + A 



(U{E,va) , 
^\V{E,va) ) W 



for weakly bound nuclei, which equally works both for 
spherical and axially deformed nuclei. In the above equa- 
tion, r and a are the coordinate of the particle in nor- 
mal and spin space, ft,(rcr, rV) and hira, r'cr') denote the 
particle-hole and the particle-particle (hole-hole) com- 
ponents of the single-particle Hamiltonian, respectively, 
U (rcr) and V{r(j) the upper and the lower components of 
the single-quasiparticle wave function, and A is the chem- 
ical potential For simplicity of notation, the isospin 
index q is omitted in Eq. ([1]), but, of course, we solve the 
CHFB equation for coupled systems of protons and neu- 
trons. In this section, we outline the calculational scheme 
and details will be presented in the succeeding sections. 

The proposed method to solve the CHFB equations, 
abbreviated HFB/PTG, consists of the following steps: 

(1) One starts with spherical or deformed HFB cal- 
culations in the HO basis (HFB/HO). This provides a 
good approximate solution for the HF potential and the 
effective mass. 

(2) One considers a HF potential and an effective 
mass for each Ij subspace, and fits the associated shifted 
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PTG potential to them when the HF potential possesses 
bound or narrow resonant states in this £j subspace (see 
Sec. IIII A[) . If no such states appear in the HF ij spec- 
trum, a set of Bessel/Coulomb wave functions [2J] is se- 
lected for the ij partial wave basis. 

(3) One diagonalizes the HFB eigenvalue equations in 
the basis composed of the PTG and Bessel/Coulomb 
wave functions. This step continues until self-consistency 
is achieved. 

The use of the Bessel/Coulomb wave functions in step 
(2) occurs for partial waves of high angular momentum, 
for which the centrifugal part becomes dominant. As 
no resonant structure can appear therein in the real HF 
continuum, Bessel/Coulomb wave functions provide a nu- 
merically stable set of states for this partial wave. For 
the generation of Coulomb wave functions, one can use 
the recently published C++ code ^ or its FORTRAN 
alternative '26] . A complete set of wave functions is thus 
formed, which will be used as a basis to expand the HFB 
quasiparticle wave functions. 

The necessary truncation of the basis in step (3) im- 
plies that spurious effects may eventually appear at very 
large distances, where both the particle density p and the 
pairing density p are very small. Consequently, quasipar- 
ticle wave functions have to be matched to their exact 
asymptotics at moderate distances as it is explained fur- 
ther in Sec. IIVI In addition, special care must be taken 
to calculate matrix elements due to the presence of non- 
integrable scattering states (see Appendix IB|) . 

When the HF mean-field resulting from the HFB/HO 
calculation in step (1) is deformed, there are several ways 
to extract the HF potential for each £j subspace to be 
used in step (2). Because it is used just as a gener- 
ator for the complete PTG basis, its choice will have 
little effect on the final HFB solution, however. In the 
present calculation, we therefore adopt a simple proce- 
dure; the particle-hole part of the HFB/HO potential 
and the HFB/HO effective mass are used in step (2) af- 
ter averaging their angular and spin degrees of freedom. 
The resulting HF potential is spherical and the same for 
all £j subspaces. In such a case, the effect of the spin- 
orbit splitting is not taken into account in the stage of 
constructing the PTG basis but it is of course taken into 
account in step (3). This implies to consider a basis gen- 
erated by a spherical potential, which might seem inefh- 
cient in the case of large deformation, for which deformed 
bases are more appropriate, as is done with the HO and 
THO bases. The deformed nuclei considered in this paper 
are nevertheless fairly reproduced within this framework 
(see Sec.|V|. If necessary, it is possible to generate a de- 
formed basis by diagonalizing the deformed HF potential 
within the PTG basis, which can then serve as a particle 
basis for the HFB problem. 



III. GENERATION OF BASIS 
A. PTG potentials fitting procedure 

The PTG potential has four parameters A, s, v and a, 
which have to be determined in each £j subspace (see 
Appendix For this purpose, we use the spherical HF 
potential and effective mass in a given £j subspace. 




FIG. 1: (color online) The shifted PTG potential, the HF po- 
tential calculated with the SLy4-force and the unshifted PTG 
potential for neutrons in ^''Ni. HF and shited PTG potentials 
to which centrifugal part is added are provided as well, and 
the energies of 0(?7/2 levels for each potential are indicated. 
All data respectively associated to HF, shifted and unshifted 
PTG potentials are respectively shown in solid, dashed and 
dotted lines. 

The PTG mass parameter a is obtained from the re- 
quirement that the PTG and the HF effective masses are 
the same at the origin. One first adds the centrifugal 
term V£(^+i) oc £{£ + l)/f^ to the nuclear plus Coulomb 
potential, Vn + Vc, and determines the height Ei, of the 
centrifugal (plus Coulomb) barrier. Then, one adds Eb to 
the PTG potential; the resulting potential may be called 
the shifted PTG potential. The parameters A and are 
fitted in such a way that the difference between the 
shifted PTG potential and the HF potential is minimal. 
Note that s is directly obtained from A and v values dur- 
ing the fit, as it is determined by way of the property that 
the PTG potential of parameters A, s, v and a for r ^ 
is equivalent to times the PTG potential of parame- 
ters A, s = 1,1/ and a. The reason why we use the barrier 
height E\i in our fitting procedure will become apparent 
by an illustrative example presented below. 

To test the fitting procedure and the quality of the re- 
sulting PTG basis we performed GHF calculations in the 
coordinate space for the spherical nucleus ^''Ni. Let us 
examine the quality of single-particle energies and wave 
functions resulting from the shifted PTG potential by 
comparing them with the GHF energies and wave func- 
tions for bound and resonance states. 
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Figure [T] illustrates the PTG fitting procedure and 
compare the results with the GHF ones taking the neu- 
tron 0^7/2 level as an example. It is seen that the energy 
of the bound 0.97/2 state in the original (unshifted) PTG 
potential (horizontal dotted line) become positive after 
being shifted with Eb (horizontal dashed line) and its po- 
sition agrees in a good approximation with the resonance 
energy obtained by the GHF calculation (horizontal solid 
line). This is due to a special feature of the PTG poten- 
tial, for which the centrifugal potential decreases expo- 
nentially and not as for r +00 (see App.|^. This 
implies that the centrifugal -I- shifted PTG potential goes 
very quickly to the constant value, Ei,, for r +00. 

In this way, the PTG treatment replaces the GHF res- 
onance with a weakly bound PTG state whose wave func- 
tion will be very similar in the nuclear region. Approx- 
imating resonant states by weakly bound states in our 
framework resembles the standard two-potential method 
described in Ref. [13]. Thus, one can expect that the 
fitted PTG potential provides a rapidly converging basis 
for solving the HFB equations. 

In fact, it is not necessary to find the PTG poten- 
tial that exactly minimize the difference with the HF 
potential. As the PTG potential is used as a basis gen- 
erator, slight differences with the exact minimum lead 
only to slightly different bases states to expand the HFB 
quasiparticle wave functions, preserving its rapidly con- 
verging properties. Thus, one can take rather large steps 
for the A, i' variations and few radii for the evaluation 
to save computer time, keeping the quality of the basis 
essentially the same. 



B. Single-particle energies 

Single-particle energies and widths for neutrons in ^^Ni 
obtained by the GHF calculations are compared with the 
PTG energies in Table HI One can clearly see the follow- 
ing facts. 

Firstly, the overall agreement between the GHF and 
the shifted PTG energies is good, which means that the 
PTG potential is flexible enough to reproduce the main 
features of the HF potential. 

Secondly, all narrow GHF resonances are represented 
as weakly bound PTG states with upward shifted PTG 
energies. This is very important because the HFB upper 
(lower) components of quasiparticle states are likely to 
have large overlaps with unoccupied (occupied) weakly 
bound and narrow resonance states. 

We note that the GHF states whose width is larger 
than 1 MeV, as a rule, are not converted to bound PTG 
states. This is not important, however, because scat- 
tering states do not exhibit rapid changes in the energy 
region of broad resonances. The broad resonance region 
can indeed be well represented in terms of the continuum 
basis states. 



TABLE L Neutron GHF levels in ^''Ni calculated with the 
SLy4 Skyrme-force and the surface-type delta pairing interac- 
tion (see Sec. V for the parameters used), which are compared 
with the PTG estimates. All energies are given in MeV while 
the width F is given in keV. 
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FIG. 2: (color onUne) The PTG (dashed lines) and GHF (solid 


lines) wave functions for various resonance states. 




As illustrated in Fig. [2] narrow GHF resonant states 


bear larg 


e overlaps with their 


associated PTG bound 


states. Hence, the GHF resonant structure present in 


the HFB quasiparticle wave 


functions will be sustained 


by the PTG bound states, thus reducing the coupling to 


the PTG scattering continuum. 






An exa 


mple indicating the quality of the bound single- 


particle wave functions resulting from the fitting PTG 
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FIG. 3: (color online) The PTG (dashed lines), GHF (solid 
lines) and HO (dotted lines) wave functions including the 
asymptotic region for the bound 0si/2,lsi/2 and 2si/2 neu- 
tron states both in normal scale (lower panel) and logarithmic 
scale (upper panel). 



procedure is shown in Fig. [3] for the bound 0si/2,lsi/2 
and 2si/2 neutron states. In this case, nuclear potential 
has no centrifugal barrier, so that the PTG and the HF 
potentials possess the same asymptotic behavior. Very 
good agreement between the PTG (dashed lines) and the 
GHF (solid lines) wave functions is thus not surprising. 
The upper panel in Fig.[3]shows the asymptotic region in 
logarithmic scale where HO wave functions (dotted lines) 
are also given as a reference. Their Gaussian asymptotics 
cannot reproduce even approximately the exponential de- 
crease of the PTG and GHF wave functions. 
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FIG. 4: (color online) The PTG (dashed lines) and GHF (solid 
lines) wave functions of the neutron continuum s-states cal- 
culated with energies of 0.118 MeV, 9.996 MeV and 66.119 
MeV. 



Neutron continuum s-states are illustrated in Fig. |31 



which are properly reproduced as well by the scattering 
states for the PTG potential. In the cases when a cen- 
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FIG. 5: (color online) The PTG (dashed lines) and the GHF 
(solid lines) wave functions for the neutron continuum ^3/2- 
states calculated at the same energies as in Fig. |4] 

trifugal (and/or Coulomb) barrier exists, as illustrated in 
Fig. [5] for (^3/2 states, different phase shifts develop in the 
PTG and GHF continuum states, as the PTG potential 
bears no barrier at large distance. 

IV. QUASIPARTICLE WAVE FUNCTIONS IN 
THE ASYMPTOTIC REGION 

The necessary truncation of the basis implies that spu- 
rious effects will eventually appear at very large radius, 
where both the particle density p and the pairing density 
p are very small. Consequently, quasiparticle wave func- 
tions have to be matched with their exact asymptotics 
at moderate distance, where the asymptotic region has 
been attained and densities are still large enough for ba- 
sis expansion to be precise. Below we explain how the 
matching procedure is done for axially deformed nuclei. 

In order to deal with the asymptotics of quasiparticle 
wave functions, wc make partial wave decomposition of 
them: 

a tj 

(2) 

a ej 

where the subscript k specifies the quasiparticle eigen- 
states together with the magnetic quantum number m 
which is always conserved for both spherical and axially 
symmetric nuclei; (r) are the PTG or Bcsscl/Coulomb 
wave functions; Uj^^^ and V^^^ are the HFB expansion co- 
efficients; ujfj^{r) and V^^\r) are the radial amplitudes 
with r — |r| for the (£j) partial wave; y^l{il.) denotes 
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a product wave function where the spherical harmonics 
with the angular variables f2 and the orbital angular mo- 
mentum i is coupled with spin to the total angular mo- 
mentum J. 

The partial wave amplitudes, ul^'^ (r) and vjlH^ (r), de- 
fined above involve summation over all quantum numbers 
except the angular momenta £ and j. In the spherical 
case, the sums reduce to a single element as i and j are 
good quantum numbers. In the asymptotic region, only 
Coulomb and centrifugal parts remain from the HFB po- 
tentials, so that one can continue the quasiparticle wave 
functions via their partial wave decompositions and de- 
cay constants A;„ and h^'- 



and discretized with Nf^j continuum states per partial 
wave in the interval [0 : k^ax]- Figure [S] shows that the 



km 



(3) 



2m 



where E denotes the quasiparticle energy, A the chem- 
ical potential, H^^^ the Hankcl (or Coulomb) functions, 



km 



km 



•q being the Sommerfeld parameter, and 

and Dj,^ are constants to be determined. Matching is 
performed using Eq. Q at a radius i?o in the asymp- 
totic region where the basis expansion is precise, so that 

^kni~^ ^ ^km~ ^'^'^ ^hm^ comc forward by continuity. 
The value of Rq is typically of the order of 10 fm. 



V. NUMERICAL EXAMPLES 

We have made a feasibihty test of the HFB/PTG 
method for spherical Ni isotopes close to the neutron 
drip line and for deformed neutron-rich nuclei ^^'^Zr and 
^°Mg. All calculations were done using the SLy4 den- 
sity functional. For the pairing interaction, we use the 
surface-type delta pairing with the strength Iq = —519.9 
MeV fm'^ for the the density-independent part and ig — 
— 37.5tg MeV fm^ for the density-dependent part with 
a sharp energy cut-off at 60 MeV in the quasiparticle 
space. They have been fitted to reproduce the neutron 
pairing gap of ^^°Sn. These values are consistent with 
those given in Ref. [32]; the slight difference is due to 
different cut-off procedures, sharp cut-off in our case and 
smooth cut-off in Ref. 32]. Below we discuss the major 
features of the result of calculation. Wc also make a de- 
tailed comparison between the HFB/PTG and HFB/Box 
calculations in the spherical case. 



A. Spherical nuclei 

Let us first examine how the result of calculation de- 
pends on the truncation of the basis. Indeed, the basis 
has to be truncated at a maximal linear momentum kmax , 




FIG. 6: (color online) Dependence on kmax of the neutron 
density p„ and the neutron pairing density p„ calculated for 
**''Ni by the HFB/PTG method. 

use of values larger than kmax — 3 fm^^ does not change 
the results. Accordingly, in calculations for spherical nu- 
clei, we use kmax = 5 fm~^ and discretize the continuum 
with Nf i = 60 scattering states per partial wave (see 
Ref. 23] for its justification). 

Results of the HFB/PTG calculation for a set of bench- 
mark Ni isotopes close to the neutron drip line are pre- 
sented in Table |TT1 Fig. [7] and Fig. [H where results of 
the HFB/Box calculation are also shown for compari- 
son. The Ni isotopes are spherical with pairing in the 
neutron channel only. We see immediately a remark- 
able agreement between the results of the HFB /PTC and 
HFB/Box calculations. The difference in total energies 
is less than 85 keV and the proton rms radii agree almost 
perfectly, while the neutron ones are slightly different by 
less than 0.003 fm. Similarly good agreement is obtained 
for all other energy counterparts. The good agreement 
in the ground state characteristics evaluated by the two 
different approaches is not surprising if one compares the 
density distributions shown in Fig. [7] and Fig.[8l In these 
figures, the neutron and proton densities, p„ and pp, and 
the neutron pairing density /5„ are plotted both in nor- 
mal (left column) and logarithmic (right column) scales. 
The agreement is almost perfect in the whole range of r 
except at the box boundary where the HFB/Box densi- 
ties vanish due to the boundary conditions (however not 
seen in Fig [8]). This agreement is striking considering the 
significant impact of the continuum for these nuclei and 
the fact that the HFB /PTG calculations are nevertheless 
performed using the basis expansion method. 

Special attention has to be paid to the agreement for 
the pairing quantities. Interestingly, the pairing gap A„ 
increases as one approaches the drip line, indicating the 
important role of the pairing correlations in the con- 
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TABLE II: Results of the HFB/PTG calculation for ground state characteristics of Ni isotopes close to the neutron drip line, 
which are compared with results of the HFB/Box calculation. The SLy4 functional and the surface-type delta pairing 20] are 
used. The rms radii are in fm and all other quantities are in MeV. Proton chemical potential Ap is not provided as pairing 
correlations vanish in the proton space. 





86 Ni 


88 Ni 


90Ni 


HFB/Box HFB/PTG 


HFB/Box HFB/PTG 


HFB/Box HFB/PTG 


HFB/Box HFB/PTG 



An 


-1.453 


-1.429 


-1.037 


-1.029 


-0.671 


-0.661 


-0.342 


-0.329 


Tn 


4.451 


4.450 


4.528 


4.526 


4.603 


4.602 


4.677 


4.674 


rp 


3.980 


3.981 


4.001 


4.001 


4.021 


4.021 


4.043 


4.043 


An 


1.481 


1.532 


1.667 


1.658 


1.790 


1.780 


1.899 


1.892 




-30.70 


-30.60 


-36.52 


-35.92 


-41.98 


-41.187 


-47.158 


-46.233 




1084.53 


1085.95 


1118.65 


1118.63 


1150.71 


1150.64 


1182.52 


1182.66 


Tp 


430.47 


430.240 


425.99 


426.01 


421.71 


421.72 


417.38 


417.37 


ipso 


-63.379 


-63.177 


-61.679 


61.707 


-59.558 


59.681 


-56.898 


-57.889 


jpCoul 


132.94 


132.90 


132.26 


132.246 


131.571 


131.578 


130.947 


130.886 


^exc 


-10.138 


10.136 


-10.084 


-10.085 


-10.033 


10.033 


-9.980 


-9.980 


Etot 


-654.89 


-654.914 


-656.933 


-656.877 


-658.167 


-658.084 


-658.665 


-658.608 



tinuum. This result is somehow different from that of 
Ref. [29! obtained by an alternative HFB calculation in 
the coordinate space for the same set of nuclei but it is in 
agreement with the estimates from Ref. [sO] • In Fig. [51 
the scaling function of the THO basis is calculated with 
the method described in Ref. [1^, for which the quasi- 
exact density provided by the HFB/PTG calculations 
is used, and 16 THO shells are taken into account for 
each partial wave. This implies virtually optimal results, 
and it has been checked that densities obtained from the 
HFB/Box and HFB/THO methods are almost identical 
up to 20 fm. On the other hand, pairing densities given 
by the THO calculations are not exactly the same with 
those of the HFB/PTG and HFB/Box calculations, as 
can be seen from Fig. [HI While pairing densities calcu- 
lated with both methods for ^"'Ni and ^"Ni are very close, 
those for *^Ni and *^Ni exhibit visible differences, espe- 
cially for ^^Ni, for which pairing energies differ by about 
4 MeV. Asymptotic properties of pairing densities calcu- 
lated with the THO basis are also not well behaved after 
15-20 fm, where they saturate instead of decreasing expo- 
nentially. This indicates that THO basis calculations are 
not always devoid of inaccuracies, even at the spherical 
HFB level. 



B. Axially deformed nuclei 

In the case of axially deformed nuclei, few HFB/Box 
calculations are available to check the HFB/PTG re- 
sults. We consider the well-deformed nucleus ^^°Zr (de- 
formation P « 0.4), already studied in Ref. [2l[ and two 
states with different deformations for the drip line nu- 
cleus ^°Mg. We use therein kmax = 4 fm~^ and Nij = 30 
for all partial waves. 

Table IIHI compares the three approaches with respect 



to ground state properties of ^^"Zr. In general they 
yield similar values. The differences seen in Table Hill are 
partially due to different structure of the model spaces 
adopted and the associated fitting of the pairing strength. 

TABLE HI: Comparison of ground state properties of ^^°Zt 
calculated with the HFB/Box, HFB/PTG and HFB/THO 
approaches. The rms radii are in fm, quadrupole moments 
are in barn, and all other quantities are in MeV. 





HFB/Box 


HFB/PTG 


HFB/THO 


Qtot 


12.088 


12.53 


12.303 


A„ 


0.480 


0.626 


0.562 


j^pair 


-1.53 


-3.015 


-2.05 




4.82 


4.836 


4.831 


rp 


4.55 


4.560 


4.556 


Etot 


-893.93 


-893.952 


-893.711 



Proton and neutron densities for nuclei ^^°Zr and ''"Mg 
are displayed in Fig. [HI with comparison with THO re- 
sults (circles) for ^^'^Zr, in normal scale (left column pan- 
els) and logarithmic scale (right column panels). Associ- 
ated pairing densities are shown in Fig. IIOI 

While agreement between the PTG and THO densities 
for -'^^"Zr is good in normal scale, we can notice discrep- 
ancies in asymptotic properties, which are visible from 
the figure in logarithmic scale (see Fig. It is obvious 
that all densities calculated with the THO basis eventu- 
ally follow the common asymptote dictated by the scaling 
function, while they are well reproduced with use of the 
PTG basis. This comparison also confirms the presence 
of deformation effects even in the far asymptotic region. 

The middle and bottom panels in Figs. [HI and [TUl il- 
lustrate the HFB/PTG normal and pairing densities for 
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FIG. 7: (color online) The neutron densities p„ and proton densities pp both in normal (left-hand side) and logarithmic (right- 
hand side) scales. Results of the HFB/Box calculation are displayed by solid lines, while those of the HFB/PTG calculations 
by open circles and dashed lines. The HFB/HO densities are also indicated by dotted lines in the right panels for comparison. 



two states with different deformations in the drip line 
nucleus ^"Mg. These states possess pairing correlations 
in both neutron and proton channels. The prolate and 
oblate states lead to asymptotic neutron densities which 
are very close, as seen from the middle and bottom right 
panels in Fig. [9l 



VI. CONCLUSIONS 

We have proposed a new method of the CHFB calcu- 
lation for spherical and axially deformed nuclei, which 
properly takes the continuum into account. The method 
combines configuration-space diagonalization of the HFB 
Hamiltonian in the complete set of analytical PTG and 
Bessel/Coulomb wave functions with a matching proce- 
dure in the coordinate space which restores the correct 
asymptotic properties of the HFB wave functions. The 
PTG potential is chosen to fit the nuclear HF potential 
and effective mass. The resulting PTG wave functions 
are close to the bound and continuum states of the related 
HF potential while the resonance states are substituted 
by the bound PTG states with shifted single-particle en- 
ergies. Partial waves of high angular momentum are very 
well represented by Bessel/Coulomb wave functions. 



The main results of the present work are twofold: 

First, we have obtained a new scheme (HFB/PTG) to 
solve the CHFB equations as a promising tool for large 
scale calculation; its performance is comparable, some- 
times even better, to that of the HFB/THO code, for 
example. It properly takes the nuclear continuum into 
account and therefore could be used for precise density 
functional calculations for nuclei close to the drip lines. 
This HFB/PTG method can also be used to provide ac- 
curate quasiparticle wave functions for microscopic cal- 
culations of dynamics beyond the nuclear mean-field ap- 
proximation, as for example, the QRPA calculations for 
deformed nuclei. 

Second, the fact that the HFB/PTG calculation re- 
produces the results of the coordinate-space HFB calcu- 
lation with the box boundary condition (HFB/Box) even 
for nuclei up to the neutron drip lines is important. This 
result indicates the validity of the HFB/Box calculation 
which is widely used, although its validity is sometimes 
questioned when it is applied to the drip-line phenomena 
where continuum effects are crucially important [29l |. 

The inclusion of resonant structure in the basis is cru- 
cial for the success of the HFB/PTG approach. Our 
test calculations indicate significant disagreement with 
the HFB/Box result if the PTG bound states represent- 
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FIG. 8: (color online) The neutron pairing densities p„ in normal (left-hand side) and logarithmic (right-hand side) scales. 
There are no pairing correlations in the proton channel. Results of the HFB/Box and HFB/PTG calculations are displayed 
both by solid lines, as they are almost indistinguishable, while HFB/THO pairing densities are represented by dashed lines. 



ing the resonant GHF states are removed from the basis: 
in their absence, the pairing densities are overestimated 
in the surface region, while particle densities are slightly 
underestimated in the inner region. This means that the 
resonance states significantly contribute to the total en- 
ergy through both the particle-hole and particle-particle 
channels. Their contributions to the pairing correlation 
energy are evaluated to be about 2-3 MeV for the case of 
Ni isotopes close to the neutron drip line. 

A more complete investigation of the importance of the 
HFB resonance states could be made by a detailed com- 
parison with the result of the exact Gamow-HFB calcu- 
lation. Such an analysis is in progress for spherical nuclei 
and will be reported in the near future [3l| . 
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APPENDIX A: PTG BASIS 

1. PTG potential 

The one-body Hamiltonian for the exactly solvable 
PTG model reads: 

H ( ^ ^ d ^ i{e+i) \ 

2mo V dr fi{r) dr r'^ii{r) J 
+ VpTG{r) (Al) 

with mg the particle free mass, r is the radial coordinate 
(in fm), /x(r) its dimensionless effective mass (the full ef- 
fective mass is mo M ('')), ^ its orbital angular momentum 
and VpTG is the PTG potential. The potential Vptg(^) 
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FIG. 9: (color online) The neutron and proton densities of the prolately deformed nucleus ^^"Zr {(3 — 0.40), respectively 
calculated by the HFB/PTG (respectively solid and dashed lines) and HFB/THO (circles) methods in normal (top left) and 
logarithmic (top right) scale. They are given along the long and short axes of deformation, easily identified from the figure. The 
neutron and proton densities of ^''Mg calculated by the HFB/PTG method for two states with different deformations (oblate 
/3 — —0.09 and prolate (3 — 0.26) in normal (middle and bottom left) and logarithmic (middle and bottom right) scale are also 
provided with the same line convention. 



and the effective mass fi{r) are written: 



/i(r) 
VpTcir) 



i-a(i-r), 



X (y^(r) + V,{r) + V,{r)), 



(A2) 



(A3) 



where s is the scaling parameter, the potential part 
issued from the effective mass, Vi its ^-dependent part 



and Vc its main central part, defined by 
V^,{r) = [l-a + [a(4-3A2) 
- 3(2-A2)]y2 

■(l-2/2)(i + (A2-l)y2) 



^(r)2' 



T4(r) = ^(£+1) 
1 



y 



r > 0, 



A2-1 



(2-(7-A2)y2 



(A4) 



(A5) 



(A6) 



- 5(A2-l)y^)]. 

The quantities VpTG{f) and ^{r) depend on an implicit 
function y = y{r) defined in the following way: 



A^s r — arctanh(y) + VA^ — 1 arctan(\/A^~~T y) 

(A7) 
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FIG. 10: (color online) Same as in Fig. |9] but for pairing densities and without HFB/THO results. Proton pairing density is 
not represented for ^^°Zr as it is negligible therein. 



so that < y < 1 for < r < oo. 

The numerical solution of Eq. (IA7p by way of New- 
ton/bisection methods is stable but one should take spe- 
cial care at large distances when y becomes closely equal 
to one. For example, this can be done by rewriting 
Eq. (|A7|1 . introducing the new variable x = arctanh(y): 



A^s r = X 



+ VA^ - 1 arctan(VA2 - ltanh(a;)), 



(A8) 



It is solved with respect to x with a fixed-point algorithm. 
In this region, 1—y^ should be calculated in terms of the 



7(1 



to avoid numeri- 



expression 1 — y = 4e" 
cal cancellations. 

One has to mention that, in the calculation of Vptg {'"') i 
Ve{r) is finite for all r > but is the difference of two 
diverging terms for r — )■ 0. Thus, to be precise in this 
region, Eq. (jA7P must be rewritten as a power series in 
y, so that the main diverging terms of Eq. (jA5[) cancel 
analytically. 

As seen from the equations above, there are four pa- 
rameters in the PTG model: the effective mass parameter 
a, the scaling parameter s, the parameter A determining 
the shape of the potential and the parameter v associated 
with the depth of the potential. They can take different 
values for different angular momenta £. We can use this 
freedom in order to approximate the nuclear potential for 



each ^j-subspace, as described in Sec. HH 



2. PTG states 

The PTG wave functions and eigen-energies are deter- 
mined by the Schrodinger equation for the Hamiltonian 



HpTG = E *fc(r) 

with energies 



2mo 



(A9) 



(AlO) 



where k stands for the complex linear momentum asso- 
ciated with E. 

For bound states, if they exist, the parameter v de- 
termines the maximal value rimax of the radial quantum 
number n = 0, 1, 2, rimax as the largest integer inferior 
to 

■ (All) 



2 V 2 



and defines the complex momentum 



knl — 



l-a 



(A12) 
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with 



Ani = 2n + £+-, 



A„i = K^\v+^] (1 - a) 



- [{l-a)K^-l]Al, 



(A13) 



(A14) 



For continuum states, k can take any real positive values 
from zero to infinity. 



3. PTG wave functions 

In order to express the PTG wave function $fe(r) = 
r \I/fe(r) in a closed analytical form, let us introduce the 
following three functions 



(A15) 



states. This defines all other quantities entering the equa- 
tions above. For both cases, the PTG wave functions can 
be written either as 



^k{r) Xk{r) fk{r) 



(A24) 



or as 



Mr) ^^fXk{r) {A+ f+{r)+A- f-(r)) . (A25) 



Equation (jA24ll is suitable for numerical work for small 
distances since ~^ when r ^ so that one is away 
from the pole of the hypergeometric function appearing 
at x~ = 1. Similarly, Eq. (|A25|) is applicable for large 
distances since — > when r — s- +oo and the pole 
2:+ = 1 of the hypergeometric function in Eqs. (jA16|) 
and (|A17p is avoided. 

In the case of bound states, the quantmn numbers {n£} 
are the principal quantum number n and the angular mo- 
mentum £. The constants A/", A~^,A^ entering Eqs. (|A24p 
and (|A25P are given by: 



and 



Xkir) 



where 



/+(r) = F + 1, x+) {x+f^ , (A16) 



f^{r)^F{^x-,^Ji+-p+l,x+){x+) (A17) 



':^il±£il=^(.-)^, (A18) 

\/x- + K^x+ 



N = 



I 2A2s/3(^ + I + /? -f 2n) 
{£+^+PA^il-a) + 2n) 



-f3 + n)T{e- 



r(n + i)r(/? + n + i)r(^-t-|)2' 



A+ = 



T{£- 



r(/i+)r(A.-) 



0, 



(A26) 



where r(z) is the Gamma function [24| . 

In the case of scattering states, the quantum numbers 
{k£} include the momentum k and the angular momen- 
tum £ while the associated constants A/", A'^, A^ read: 



_ 1 _ (A2 + 1)^2 
l + (A2-l)y2'^ 



i_^,^+^i±^, (A19) 



^__£+'^+P + D __£+l + [3-v 



. £+l-p + i? 
" =^-2 ' ^' 



(A20) 



, (A21) 



(A22) 



AA = 



/ r(^+)r(z.-)r(M+)r(/i-) 
27: r(/3)r(-;3)r(^ + f)2 



Ti£+ln-f3) 

r(M+)r(M-) 
T{£ + I)r{p) 

r(^+)r(^-) • 



A- = 



(A27) 



The normalization constant J\f is determined from the 
normalization condition 



<i>„iir)<i>n'i{r)dr = d„ 



(A28) 



for bound states and from the Dirac delta function nor- 
malization for scattering states: 



(i. + l/2)2-|-^2(i_ A2(l-a)), (A23) 



and F(a, b, c, z) is the Gauss hypergeometric function 

In the case of bound states, fc„; determines the mo- 
menta k which are pure imaginary (see Eq. (jA12p ). while 
they are real positive numbers in the case of scattering 



^ki{r)<^k'i{r)dr = S{k--k') 



(A29) 



All bound and scattering wave functions arc orthogo- 
nal to each other 



/"OC 

/ ^k{r)<Pk'{r)dr = 0, k ^ k' 

JO 



(A30) 
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and they form a complete basis 

71 1 

+ f 'fu{r)<fki{r')dk^6{r~r'). 
I 

One can check that at large distances 

a;^ -l + 2e-2A''s('-n)^ r +00, 

where 

A 



(A31) 



A^s ri = - larctan(VA2 - l)-log ( — 



(A32) 



(A33) 



Substituting this into Eq. (jA25[) one obtains the asymp- 
totic form of the PTG wave functions 



$fc(r) ^ c+ e^"" + c- 



(A34) 



where C'+ = TVA+g-^i and C = MA-e^^''\ defined 
by Eqs. ([JMjl and (|M7)) . 

The PTG wave functions are numerically stable and 
accurate when using Eq. (jA24[) up to ?/ < 0.99 then 
applying the form (|A25p . They accurately land onto 
their asymptotic representation of Eq. (IA34p at large dis- 
tances. 



APPENDIX B: MATRIX ELEMENTS 

Let us deal with numerical integration in r and h space. 
The integration in the r space is performed in terms of 
Nr Gauss-Legendre integration points Xi and weights Wi 
within the interval [0,i?maa:], 



0(r) $fe(r) $fc/(r) dr 







^ 0(ri) ^k{r^) ^k'(ri)w^, 



(Bl) 



where 0(r) is an arbitrary function of r and Rmax is a 
point where nuclear potential disappears. Usually a value 
Rmax = 15 fm is used. In the same way, integration 
in the k space is done in terms of Nk Gauss-Legendre 
integration points ki and weights Wk^ within the interval 
[0, k^dx] , 



0{k) $fc(r) $fe(r') dk 







potential or explicitly depends on nuclear densities or 
currents, one can safely integrate the matrix elements to 
some large but finite distance Rmax- Beyond Rmax, the 
contribution of the integral becomes negligible due to the 
presence of the densities or currents. However, it is not 
the case for the kinetic -t- Coulomb part of the Hamilto- 
nian. This operator is infinite-ranged and induces Dirac 
delta functions in the matrix elements, which have to be 
regularized directly. For this, one separates the matrix el- 
ement in two integrals, defined on the intervals [0 : Rmax] 
and [Rmax '■ +00 [. The first part is finite and treated with 
standard methods. For the second part, if one deals with 
Bessel/Coulomb wave functions, one can assume that the 
nuclear part is negligible after Rmax so that they are so- 
lutions of the asymptotic HE equations. Hence, one ob- 
tains: 



+00 



Ua{r)h{r)u/3{r) dr 



klySai3-J Ua(r)u^(r) drj (bound) 

kl is{ka - kfj) - J Ua{r)u/3{r) drj (scat) 



~kt I Ua{r)ufj{r) dr (mixed) 
'0 



(B3) 



0{h) ^kAr) ^kAr') u-fe., 



where h{r) is the HE potential which reduces to the 
kinetic + Coulomb Hamiltonian asymptotically. Here, 
"bound" ("scat") means that both a and (3 states are 
bound (scattering) and "mixed" means that a is bound 
and 13 scattering or vice-versa. The Dirac delta with 
a discretized basis becomes Sap/wk^ with Wk^ being 
the Gauss-Legendre weight associated to the discretized 
value fca, so that its implementation is immediate; since 
all integrals are finite, they pose no problem. When the 
PTG basis states are used instead of the Bessel/Coulomb 
wave functions, it turned out that it is numerically precise 
to disregard the Coulomb/centrifugal part of the Hamil- 
tonian after Rmax, so that Eq. (jB3p is the same for both 
the PTG and Bessel/Coulomb wave functions. Indeed, 
Eqs. (|M2)) and (|X34)) imply that the PTG wave func- 
tions behave asymptotically like neutron waves functions 
of angular momentum t — Q. The above seemingly crude 
approximation can, in fact, be mathematically justified. 
The HEB matrix evaluated using such a procedure con- 
verges weakly to the exact HEB matrix for Rmax +00 
[3^ . This means that the HEB matrix elements depend 
on Rmax asymptotically, some of them even diverging 
with Rmax — *■ +00, whereas its eigenvalues and eigenvec- 
tors converge to a finite value. 



(B2) 



where 0{k) is an arbitrary function of k. 

Radial integrals must be calculated cautiously due to 
the presence of non-integrable scattering states in the 
basis. When the radial operator represents the nuclear 
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